Densification effects on the Boson peak in vitreous silica: 
a molecular-dynamics study 
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We perform classical molecular-dynamics simulation to 
study the effect of densification on the vibrational spectrum 
of a model silica glass. We concentrate this study on the so- 
called Boson peak and compare our results, obtained from a 
direct diagonalization of the dynamical matrix, with experi- 
mental Raman data. We show that, upon densification, the 
position of the Boson peak shifts towards higher frequencies 
while its magnitude decreases which is in agreement with a 
recent experimental study. 

PACS numbers: 61.43.Fs, 63.50.-tx, 02.70.Ns 



I. INTRODUCTION 

Vitreous silica, the prototype material for the study of 
network forming glasses, displays a number of anoma- 
lous behaviors, both in the structural and vibrational 
characteristicall. This is also true when this material is 
exposed to high pressures. Experiments have shown that 
at room temperature glass samples compressed beyond 
« 12jpPa exhibit a permanent density increase of abc 
20 %□. Both experimentsu and numerical simulations 
have shown that this irreversible densification is mainly 
due to an increase of the average coordination number 
around the silicon atoms. An other topic concerns the 
influence of the pressure on the vibrational spectrum of 
vitreous silica. Very recently a substantial amount of ex- 
perimental work has been dedicated to the study of this 
influence and more precisely to the effect-of the densifica- 
tion on the so-called Boson peak (BP)Qij. The BP refers 
to an excess in the vibrational density of states with re- 
spect to the Debye distribution and is located generally 
around 1.5 THz. It is not clear yet if the BPrOnly shifts 
towards higher frequencies upon densificatioi£l or if this 
shift is a result of the suppression of the BPEI. With the 
use of numerical simulations using a relatively realistic 
interaction potential one should be able to give a clear 
answer and connect this answer to the structural traps.- 
formations observed in the afore mentioned studiesoQ. 
As far as we know no simulations were done on the vari- 
ation of the BP with the density. Here we present such a 
study in a model silica glass using the same pair-potential 
than Tse et alB, the "BKS" potential^. This potential 
has also been used recoatly to study pressure induced 
amorphization of quartztil and seems appropriate to de- 
scribe silica at high pressure. Our results showing a shift 
of the BP as the density increases between 2.2 and 2.67 



g/cm^ are in relatively good agreement with the experi- 
mental Raman data and indicate a disappearance of the 
BP upon deuaification. In parallel we find similar results 
as Tse et al.a concerning the structural changes around 
the silicon atoms which are therefore not shown again 
here. 



II. SIMULATIONS 

We perform classical (constant energy, constant vol- 
ume) molecular-dynamics simulations on an assembly of 
648 particles (216 Si02 molecules) packed in a cubic box 
with periodic boundary conditions. The size of the box is 
fixed in order to study samples at the following densities: 
2.17, 2.32, 2.40 and^.67 g/cm^. The particles interact via 
the BKS potentiaJl3 with the same parameters as in our 
previous studiesO. The low temperature glassy samples 
at a given density are obtained after a quench from the 
liquid state (T « 7000K) at a constant quenching rate of 
2.3 X lO^^K/s. It is worth noticing that this procedure is 
different from the one used intl where a compression at 
room temperature was used in order to obtain the high 
density samples. After the quench, the samples at zero 
temperature are relaxed during 120000 timesteps (84 ps) 
(during which the temperature rises only slightly) and fi- 
nally the vibrational spectrum g{i') is obtained by direct 
diagonalization of the dynamical matrix. This direct di- 
agonalization is the most computer time consuming since 
one has to deal with a 1944 x 1944 matrix with no a pri- 
ori symmetry. Moreover since the Coulomb interactions 
have not been cut off, this matrix can not be considered 
as sparse (to reduce the computer time it should also be 
noted that we did not calculate the eigenvectors in this 
simulation). The eigenvalues of the dynamical matrix 
can also be obtained from the Fourier transform of the 
velocity-velocity autocorrelation function but in the low 
frequency region of the spectrum (which is the one we 
are mostly interested in) this technique is generally less 
reliable. For each density two independent liquid sam- 
ples were used in order to improve slightly the statistics 
of the results. Because of the finite size of the simulation 
box we have discrete values of the frequency v which are 
located between « 0.8 and 40 THzJn agreement with a 
previous study on the same systemllj. Nevertheless since 
we are interested in the low frequency part of the spec- 
trum we will concentrate here on the frequency region 
below lOTHz. 
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III. RESULTS 

A standard way of extracting the BP from the vibra- 
tional spectrum is to plot g{v)/v^ since in the Debye 
approximation g{v) k, at low frequency. In fig.l we 
show the averaged (over 2 samples) spectrum at differeruL 
densities together with the neutron diffraction resultsB 
when available (the curves have been arbitrarily shifted 
for clarity). The solid lines in the figure are the fits of 
g{v)/v'^ by_the following "generalized Lorentzian" func- 
tional fornJlil: 



where /o, vq, m and n are adjustable parameters (this 
functional form is useful to describe the BP even though 
it has no direct physical meaning). In this figure the 
excess of modes between 1 and 2 THz is clearly visible. 
Moreover a broadening of the peak when increasing the 
density can be seen in agreement with the experimental 
results. The increase of the peak position with increasing 
density is less pronounced than in experiment because of 
finite size effects. Indeed due to the limited size of our 
simulation box the number of frequencies "available" on 
the low frequency side of the spectrum is small. This 
explains the overestimated peak position at 2.2 g/cm^ 
compared, to the experimental data as has been shown 
recentlyllj. We did the calculation for one sample at 
2.2 g/cm'^ containing 1536 particles and obtained very 
similar results except a small shift of the peak position 
towards lower frequencies. Increasing the number of par- 
ticles by a factor 2 improves barely the description of the 
peak, since the lower cut-off frequency is proportional to 
the inverse of the box size. In any case increasing the 
system size leads to a decrease ©tthe peak position and 
an increase of its intensity (seell3) while increasing the 
density leads to an increase of the peak position and a 
decrease of its intensity as can be seen in fig. 1. The finite 
size effects are certainly responsible for the differences be- 
tween the simulation and the experimental data but do 
not affect the general behavior with density. Since only 2 
densities have been considered in the neutron diffraction 
study, we decided to compare our results to the more nu- 
merous Raman data availableBEl. When doing this one 
has to be careful in order to compare things that are 
comparable. The reduced Raman intensity is related 
to the vibrational density g{h') via a frequencjj-dependent 
function C{v). Previous experimental studiesEj and a rje-, 
cent study combining experiments and MD simulationsEZI 
have indicated that G{v) « v in the region of the BP 
which implies that Ir oc g{v)/v. Hence in the following 
we use the function g{v)/v and compare its character- 
istics to the available Raman data. To put this com- 
parison on more quantitative grounds we fit the curves 
g{v)/v = /{v) by a functional form similar to the one 
given in Eq. 1, and we extract the position of the maxi- 
mum Vrriax ^s well as the intensity of the maximum Imax 



from the fits. In fig. 2a and 2b the values of Umax and 
Imax are plotted as a function of the density respectively 
and compared to experimental Raman values. In fig. 2a, 
the 2 points reproduced fromQ (open squares) which is a 
study of the vibrational properties of silica as a function 
of pressure were obtained assuming a 15 % densification 
of the high pressure sample as stated by the authors. 

IV. DISCUSSION 

The results reported in fig. 2 show clearly that the BP is 
shifted towards higher frequencies and its intensity van- 
ishes when the density of our model silica glass increases. 
In fact, compared to fig.l, not surprisingly, the relative 
variations of the peak position and intensity are more 
important when g{v)/v is used. In fig. 2a we see that the 
overall trend of the BP position as a function of density 
is correct even though our MD simulations overestimate 
this position compared to the experimental data which is 
probably due to finite size effects as stated earlier. Nev- 
ertheless it appears that our results come closer to the 
experimental data at high density even though the posi- 
tive curvature seems to be absent in the simulation (the 
error bars in fig. 2a and 2b result from the dispersion over 
the 2 samples studied here and therefore should be con- 
sidered as a rough estimate of the error). This shows 
that the shift of the BP towards higher frequencies is 
more important than the artificial drift due to the fact 
that with increasing density (decreasing system size) the 
lower cut-off frequency increases in our simulations. This 
increase would imply a worse description of the BP re- 
gion but since the BP shifts towards higher frequencies 
our description becomes (relatively) better at high den- 
sity. Another point concerns the infiuence of the rapid 
quenching rate: indeed VoUmayr et al. have shown that 
the height and to a lessepdegree the shape of the BP de- 
pend on the cooling rateEJ. Nevertheless, similarly to the 
finite size effects, these authors show that a smaller cool- 
ing rate would lead to a behavior of the BP opposite to 
the one observed with increasing density. Therefore since 
all the samples have been quenched at the same rate the 
evolution of the BP with the density can not be linked to 
cooling rate effects. In fact all the previous studies show 
that using larger system sizes and smaller cooling rates 
would lead to a better agreement of the curves plotted in 
fig. 2a especially at low density but would not affect the 
behavior of the BP with the density which is the aim of 
this study. 

Differences can also be seen between the experimental 
studies. The principal distinction between the two exper- 
imental procedures is the temperature since Sugai et al. 
performed the densification at room temperature while 
Inamura et al. compressed the samples at 700 C. This 
might explain why our low temperature simulations are 
closer to the experimental results obtained at room tem- 
perature. In order to check the effect of the temperature. 
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we performed some calculations on samples quenched at 
300K but no major differences were noticeable. Also we 
tried a "cold compression" at BOOK followed by a long 
relaxation run (200 ps) but again the spectra were sim- 
ilar to the ones reported in fig.l. The results in fig. 2b 
confirm that with increasing density the BP tends to dis- 
appear. Again we do not find the almost linear decrease 
observed in experiment at low density but such a behav- 
ior can be observed for the higher densities (this shows 
again that our description becomes better at high den- 
sities). Since we fixed the calculated and experimental 
intensity of the BP to 1 at the usual density (2.2 g/cm"^) 
for which we know that the calculated intensity is un- 
derestimated because of the finite size of the simulation 
box, we have a systematic overestimation of the relative 
intensity of the BP at higher density. Nevertheless the 
overall trend is dear and it is coherent with the data 
of Inamura et alB who recently showed not only a shift 
towards high frequencies but also a disappearance of the 
BP. In their publication these authors attributed the sup- 
pression of the low-energy dyna|mics to the shrinkage of 
the void space by densificationt^l. This analysis is co- 
herent with the soft potential model (the soft modes in 
the void space have been suppressed by shrinkage) but 
seems to contradict the analysis of the shift of the BP be- 
ing due to phonon scattering by local density|-fluctuations 
which would result in a shift of the BP onlyQ. Neverthe- 
less the microscopic model of the BP is not settled yet 
since the structural entities involved in the soft potential 
model have not been clearly identified. The connection 
between this study and the previous numerical studies 
investigating the structural changes under densification 
should permit to progress in such an identification. 

V. CONCLUSION 

We performed classical molecular-dynamics simula- 
tions combined with the diagonalization of the dynamical 
matrix to investigate the influence of the density on the 
Boson Peak in a model silica glass. Even though we used 
a relatively small system size and a rather fast cooling 
rate, the results compare relatively well with experimen- 
tal Raman data and show that upon densification the 
Boson Peak shifts to higher frequencies while its inten- 
sity decreases. These results confirm the disappearance 
of the BP in the high, density samples as shown in a recent 
experimental studja which could be due to a shrinkage 
of the void spacdlB. The quantitative agreement between 
our simulations and experiment is not perfect (it should 
be noted that differences exist also between the experi- 
mental data) but we know the reasons of this discrepancy. 
Nevertheless the BKS potential is able to give the correct 
behavior of the vibrational spectrum under compression 
and this shows once more the good quality of this pair- 
wise interaction. Of course since we did not calculate 
the eigenvectors we have no informations on the nature 



of the modes dissapearing under compression. But the 
disappearance of the excess of low-energy modes under 
densification is the first step in the clear identification of 
what the vibrating structural entities are. To perform 
this identification in a simulation one needs tools more 
accurate than the radial pair distribution function or the 
Voronoi' tessellation. Also if the connection between the 
BP and the plateau in the thermal conductivity exists 
one should see a variation of the plateau with the densitsL. 
this variation has already been reported experimentallyEj 
but not yet confirmed in a numerical simulation. 
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FIG. 1. Plot of g{v)/v'^ versus u as a. function of density: 
• : 2.17 g/cm^; o: 2.32 g/cm^ ■: 2.40 g/cm^; □: 2.67 grjcm^ 
Comparison with experimental neutron diffraction datafl (the 
intensity at 2.2 g/cm"^ has been adjusted to the corresponding 
MD intensity ): *: 2.2 g/cm^; x: 2.63 g/cm^. 
The solid lines correspond to the fit of the MD results with 
the "generalized Lorentzian" given in Eq. 1 (see text) 
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FIG. 2. (a) Variation of the position of the BP as a func- 
tion of density; (b) Variation of the intensity of the BP as a 
function of density. 
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